Weakly coupled, antiparallel, totally asymmetric simple exclusion processes 



O 

o 

(N 

wo: 
< 



43 

o 



Robert JuhasJB 

Research Institute for Solid State Physics and Optics, H-1525 Budapest, P.O.Box 49, Hungary 

(Dated: February 6, 2008) 

We study a system composed of two parallel totally asymmetric simple exclusion processes with 
open boundaries, where the particles move in the two lanes in opposite directions and are allowed to 
jump to the other lane with rates inversely proportional to the length of the system. Stationary den- 
sity profiles are determined and the phase diagram of the model is constructed in the hydrodynamic 
limit, by solving the differential equations describing the steady state of the system, analytically 
for vanishing total current and numerically for nonzero total current. The system possesses phases 
with a localized shock in the density profile in one of the lanes, similarly to exclusion processes en- 
dowed with nonconserving kinetics in the bulk. Besides, the system undergoes a discontinuous phase 
transition, where coherently moving delocalized shocks emerge in both lanes and the fluctuation of 
the global density is described by an unbiased random walk. This phenomenon is analogous to the 
phase coexistence observed at the coexistence line of the totally asymmetric simple exclusion pro- 
cess, however, as a consequence of the interaction between lanes, the density profiles are deformed 
and in the case of asymmetric lane change, the motion of the shocks is confined to a limited domain. 
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I. INTRODUCTION 



The investigation of interacting stochastic driven diffu- 
sive systems plays an important role in the understand- 
ing of nonequilibrium steady states [H, E|- As opposed 
to equilibrium statistical mechanics, phase transitions 
may occur in these systems even in one spatial dimen- 
sion Q. The paradigmatic model of driven lattice gases 
is the one-dimensional totally asymmetric simple exclu- 
sion process (TASEP) @, [1, which exhibits boundary 
induced phase transitions Q and the steady state of 
which is exactly known [7, 8]. Beside theoretical interest, 
this model and its numerous variants have found a wide 
range of applications, such as the description of vehicular 
traffic [9[ or modeling of transport processes in biologi- 
cal systems Inspired by the traffic of cytoskeletal 
motors such models were introduced where a to- 
tally asymmetric exclusion process is coupled to a finite 
compartment where the motion of particles is diffusive 
0, EE 13 • Recently, the attention has turned to exclu- 
sion processes endowed with various types of reactions 
which violate the conservation of particles in the bulk 

El EE Ei M EE EE EH, El El, EE EE EE El- 

The simplest one among these models is the TASEP with 
"Langmuir kinetics" , where particles are created and an- 
nihilated also at the bulk sites of the system [l(| — a pro- 
cess, which may serve as a simplified model for the co- 
operative motion of molecular motors along a filament 
from which motors can detach and attach to it again. For 
these types of systems, the time scale of nonconserving 
processes compared to that of directed motion and the 
processes at the boundaries is crucial. If the nonconserv- 
ing reactions occur with rates of larger order than the 
inverse of the system size L, then in the large L limit, 
they dominate the stationary state. On the contrary, 



'Electronic address: juhasz@szfki.hu 



when they are of smaller order than 0(1/L), they are 
irrelevant and the stationary state is identical to that of 
the underlying driven diffusive system. However, in the 
marginal case when the rates of nonconserving processes 
are of order 0(1/L), the interplay between them and the 
boundary processes may result in intriguing phenomena, 
such as crgodicity breaking [IE EH or the appearance 
of a localized shock in the density profile [l^ . which is 
in contrast to the delocalized shock dynamics at the co- 
existence line of the TASEP @, El- The formation of 
domain walls can be observed also experimentally in the 
transport of kinesin motors in accordance with theoreti- 
cal predictions [IE EH EH • 

Other systems which have an intermediate complexity 
compared to exclusion processes with bulk reactions and 
those coupled to a compartment are the two-channel or 
multichannel systems. In these models, particles are ei- 
ther conserved by the dynamics in each lane and interac- 
tion is realized by the dependence of the hop rates on the 
configuration of the parallel lanes l33l [343. 13a. l36l. l37l El] , 
or particles can jump between lanes [39|, ED], Q . We 
study in this work a two-lane exclusion process where 
particles move in the two lanes in opposite directions. 
Particles are allowed to change lanes and we restrict our- 
selves to the case of weak lane change rates, i.e. they 
are inversely proportional to the system size. This means 
that the probability that a marked particle changes lanes 
during the time it resides in the system is 0(1). If par- 
ticles in one of the lanes are regarded as holes, and vice 
versa, this model can also be interpreted as a two-channel 
driven system where particles move in the same direction 
in the channels and are created and annihilated in pairs. 
In the hydrodynamic limit of the model, we shall con- 
struct the steady-state phase diagram by means of ana- 
lyzing the differential equations describing the model on 
the macroscopic scale. At the coexistence line, where co- 
herently moving delocalized shocks develop in both lanes, 
which is reminiscent of the delocalized shock dynamics at 
the coexistence line of the TASEP, the density profiles are 
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studied in the framework of a phenomenological domain 
wall picture based on the hydrodynamic description. Re- 
cently, a two-lane exclusion process has been investigated 
with weak, symmetric lane change, where particles move 
in the lanes in the same direction [42j . In this model, the 
formation of delocalized shocks in both lanes has been 
found, as well. In our model, even the case of asymmetric 
lane change can be treated analytically in the hydrody- 
namic limit if the total current is zero, which holds also 
at the coexistence line. 

The paper is organized as follows. In Sec. [Til the 
model is introduced and the hydrodynamic description 
is discussed. In Sec. IIII1 the case of symmetric lane 
change is investigated, while Sec. IIVI is devoted to the 
asymmetric case. The results are discussed in Sec. Ivland 
some of the calculations are presented in two Appendixes. 



II. DESCRIPTION OF THE MODEL 

The model we focus on consists of two parallel one- 
dimensional lattices with L sites, denoted by A and B, 
the sites of which are either empty or occupied by a parti- 
cle. The state of the system is specified by the set of occu- 
pation numbers n i ' B which are zero (one) for empty (oc- 
cupied) sites. We consider in this system a continuous- 
time stochastic process where the occupations of pairs of 
adjacent sites change independently and randomly after 
an exponentially distributed waiting time. The possible 
transitions and the corresponding rates, i.e. the inverses 
of the mean waiting times, are the following (Fig. [JJ). On 
chain A, particles attempt to jump to the adjacent site 
on their right-hand side, whereas on chain B to the adja- 
cent site on their left-hand side, with a rate which is set 
to unity, and attempts are successful when the target site 
is empty. On the first site of chain A and on the £th site 
of chain B particles are injected with rate a, provided 
these sites are empty, whereas on the Lth site of chain 
A and on the first site of chain B they are removed with 
rate f3. So the system may be regarded to be in contact 
with virtual particle reservoirs with densities a and 1 — (3 
at the entrance- and exit sites, respectively. The process 
described so far is composed of two independent totally 
asymmetric simple exclusion processes. The interaction 
between them is realized by allowing a particle residing 
at site i of chain A(B) to hop to site i of chain B{A) with 
rate ioa (j^b), provided the target site is empty. 

As in the case of the TASEP with Langmuir kinetics, 
one must distinguish here between three cases, concern- 
ing the order of magnitude of the lane change rates in 
the large L limit. If the rates to a and lob are of larger 
order than 0(1/ L), then in the limit L — > oo, the inter- 
chain processes are dominant compared to the effects of 
the boundary reservoirs and the horizontal motion of the 
particles. The densities p and tt in lane A and B, respec- 
tively, are expected to be constant far from the bound- 
aries and to fulfill the relation loap(^ — tt) = ll>b^(1 — p), 
which is forced by the lane change kinetics. When the in- 
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FIG. 1: Transitions and the corresponding rates in the model 
under study. 



terchain hop rates are smaller than 0(1/L), then (apart 
from some possible special parameter combinations [201] ) 
they are irrelevant in the L — > oo limit and the sta- 
tionary state is that of two independent exclusion pro- 
cesses. An interesting situation arises if the rates u>a and 
u>b are proportional to 1/L. In this case the effects of 
boundary reservoirs and those of lane change kinetics are 
comparable and the competition between them results in 
nontrivial density profiles. We focus on this case in the 
present work, and parametrize the lane change rates as 
to A = &a/L and = Qb/L, with the constants £Ia 
and Qb- Setting the lattice constant a to a — 1/L and 
rescaling the time t as r = t/L, we are interested in the 
properties of the system in the (continuum) limit L — ► oo, 
where the state of the system is characterized by the local 
densities p(x,r) and n(x, r) on chain A and B, respec- 
tively, which are functions of the continuous space vari- 
able x £ [0, 1] and time r. Turning our attention to the 
subsystem containing lane A(B) alone, we see that the 
interchain hoppings can be interpreted as bulk noncon- 
serving processes for the TASEP in lane A(B). The bulk 
reservoir which the TASEP is connected to is, however, 
not homogeneous but it is characterized by the position 
and time-dependent density tt(x,t)(p(x,t)). Generally, 
driven diffusive systems which are combined with a weak 
(i.e. 0(1/ L)) bulk nonconserving process are described 
on the macroscopic scale specified above by the partial 
differential equation 



dp(x,r) dJ(p(x,r)) 



dr 



dx 



S(p(x,r)), 



(1) 



where S(p(x,t)) is the source term related to the non- 
conserving process and J(p) is the current as a function 
of the density in the steady state of the correspond- 
ing translation invariant infinite system without non- 
conserving processes (i.e. S(p(x,r)) = 0) Under 
these circumstances, the TASEP has a product measure 
stationary state and the current-density relationship is 
simply J(p) — p(l — p), hence the currents in the two 
lanes are given as Ja(p(x,t)) = p(x, r)[l — p(x,r)] and 
Jb(tt(x,t)) = —tt(x,t)[1 — tt(x,t)] in terms of the lo- 
cal densities. For the source terms in the two lanes, we 
may write Sa(p(x,t),tt(x,t)) = -S b (p(x,t),it(x,t)) = 
J1_b[1 — p(x, t)]tt(x, t) — VIap{x, t)[1 — ir(x, t)] since the 
lane change events at a given site are infinitely rare in 



3 



the limit L — > oo. Setting these expressions into eq. (p} 
we obtain that in the steady state, where d T p(x, r) = 
d T Tr(x,T) = 0, the density profiles p(x) and ir(x) satisfy 
the coupled differential equations 

(2 P - i)d xP + n B (i - P )n - n AP (i - tt) = o, 

(2tt - i)a K 7r + n B (i - p)7r - n A p(i -n) = o. (2) 

We mention that one arrives at the same differential 
equations when in the master equation of the process the 
expectation values of pairs of occupation numbers {riirij} 
are replaced by the products (rii)(rij), and afterwards it 
is turned to a continuum description with retaining only 
the first derivatives of the densities and neglecting the 
higher derivatives which are at most of the order 0(1/ L) 
almost everywhere. 

For the stationary density profiles the boundary con- 
ditions p(0) — 7r(l) = a and p(l) = n(0) = 1 — (3 are 
imposed. In fact, we shall keep these boundary condi- 
tions only for a, /3 < 1/2; otherwise, we modify them 
for practical purposes at the level of the hydrodynamic 
description. The reason for this is the following. In the 
domain a, (3 > 1/2 of the TASEP, the so-called maximum 
current phase, the current is limited by the maximal car- 
rying capacity in the bulk, J =1/4, which is realized at 
the bulk density p = 1/2 @, Q ■ In this phase, boundary 
layers form in the stationary density profile at both ends, 
where the density drops to the bulk value p — 1/2. Sim- 
ilarly, in the case of the TASEP with Langmuir kinetics, 
if the entrance rate a exceeds the value 1/2, then in the 
density profile dictated by the reservoir at the entrance 
site, a boundary layer develops, where the density drops 
to 1/2. The width of the boundary layer is growing sub- 
linearly with L (2p| . such that in the hydrodynamic limit, 
it shrinks to x = and lirn^o p(x) = 1/2 holds, inde- 
pendently of a, which influences only the shape of the 
microscopic boundary layer. These considerations apply 
also to the present model at both ends and for both lanes. 
Therefore, in order to simplify the treatment of the prob- 
lem at the level of the hydrodynamic description, we use 
the effective boundary conditions 

p(0) = tt(1) = a, p(l) = tt(0) = 1 - b, (3) 

where a = min{a, 1/2} and b = min{/3, 1/2}. However, 
we stress that, although, the profile propagating from 
e.g. the left-hand boundary/?; (x) is continuous at x = 
according to the effective boundary conditions ([3]) for 
a > 0, a boundary layer forms on the microscopic scale. 

In addition to the boundary layers related to the maxi- 
mal carrying capacity in the bulk, the stationary density 
profiles may in general contain another type of bound- 
ary layer of finite width or a localized shock in the bulk, 
where the density has a finite variation within a region 
the width of which is growing sublinearly with L [T^. Il8| . 
This leads to the appearance of discontinuities in p(x) 
and it(x) in the hydrodynamic limit, either in the bulk 
< x < 1 in the case of a shock or at x = 0, 1 in the case 
of a boundary layer. This is in accordance with the fact 



that, in general, there does not exist a continuous solu- 
tion to the two first order differential equations, which 
fulfills all four boundary conditions. Apart from some 
special parameter combinations, there is one discontinu- 
ity in each lane, which is either in the bulk (a shock) or at 
x = 0, 1 (a boundary layer). The location of the disconti- 
nuity is determined by the requirement that the currents 
in both lanes Ja(p(x)) and Jb(p(x)) must be continuous 
functions of x in the bulk < x < 1 [lj], [H[ • This follows 
from that the width of the shock region is proportional 
to -\/L, thus the rate of a lane change event is vanishing 
there in the limit L — > oo. This condition permits only 
such a shock which separates complementary densities on 
its two sides, i.e. p and 1 — p in lane A or tt and 1 — 7r 
in lane B. The position of the shock x s e.g. in lane A is 
thus given implicitly by the equation pi(x s ) = 1 — p r (x s ), 
where pi[x) and p r (x) are the solutions on the two sides 
of the shock. For the detailed rules on the stability of 
the discontinuity at x = 0, 1 see Ref. fl~7j j . 

Subtracting the two differential equations yields the 
obvious result that the total current 

J = p(s)[l-/j(a;)]-7r(aj)[l-7r(aj)] (4) 

is a (position independent) constant. This relation makes 
it possible to eliminate one of the functions, say, tv(x) 
and to reduce the problem to the integration of a single 
differential equation 



dp P -[\±^{p-\Y + J][K{l-p)+p] 

dx~~ nA 2p~l ' (5) 

where we have introduced the ratio of lane change rates 
K = Hb/ and the signs in front of the square root are 
related to the two solutions tt+{x) > 1/2 and 7r_(x) < 
1/2 of the quadratic equation Disregarding the sim- 
ple case K = 1, there are two difficulties about this equa- 
tion. First, the solution depends on the current J as a pa- 
rameter, which itself depends on the density profiles and 
is a priori not known. Fortunately, apart from two phases 
in the phase diagram, p(x) and tt(x) simultaneously fit to 
the boundary conditions either at x = or x = 1, con- 
sequently, the current is exclusively determined by the 
entrance- and exit rates as J = a(l — a) — b(l — b). In 
the remaining two phases, the functions p(x) and n(x) 
meet the boundary conditions at the opposite ends of 
the system. Here, one may solve eq. ([5]) iteratively un- 
til self-consistency is attained. Second, even in the case 
when J is known, eq. ([5]) cannot be analytically inte- 
grated in general, except for the case when the current is 
zero. This is realized in three cases, two of which are re- 
lated to the symmetries of the system. We discuss these 
possibilities in the rest of the section. 

As the two entrance- and exit rates were chosen to be 
identical, the obvious relation holds when the rates £Ia 
and Qb are interchanged: 

p(x; a, (3, Q a ,^b) = tt(1 - x; a, (3, Cl B ,VL A ), (6) 
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where the dependence of the profiles on the four parame- 
ters a, f3,Q, a and Ob is explicitly indicated. This relation, 
together with eq. ((4|) implies that the current changes 
sign if Oa and Ob are interchanged. Thus Oa = Ob 
implies J = 0, that holds apparently since none of the 
chains is singled out in this case. 

As a consequence of the particle-hole symmetry of the 
model, we have the relation 

p(x; a, /3, Oa, O b ) = 1 - n(x; f3, a, A ,0 B ). (7) 

Using eq. (01, it follows that the current changes sign 
when a and (3 are interchanged, so it must be zero for 
a — j3. Alternatively, this can be seen by interchanging 
particles and holes in one of the chains, which results in a 
two-channel system where particles move in the channels 
in the same direction and particles are created and anni- 
hilated in pairs at neighboring sites of the two chains with 
rates Oa and Ob, respectively. Particles are injected and 
removed in the channels with the same rate, hence the 
channels arc equivalent. Since the currents of particles 
and holes are equal, the total current must be zero. 

The third parameter regime where the current is zero 
is the domain a, (3 > 1/2. Here, as aforesaid, the density 
profiles and the current are independent of a and (3 in the 
hydrodynamic limit. Since the current is zero for a = (3 
it follows that J = in the whole domain a, (3 > 1/2. 

III. SYMMETRIC LANE CHANGE 

We start the investigation of the model with the simple 
case of equal lane change rates (Oa = Ob = O), where 
the solutions to the hydrodynamic equations are analyt- 
ically found and some general features of the model can 
be understood. Since the total current is zero, either 
p(x) = tt(x) or p(x) = 1 — Tt(x) must hold. Substituting 
the former relation into eq. @ yields 

Pe{x) = n e (x) = const, (8) 

whereas the latter gives 

Pc(x) = 1 — k c (x) — Ox + const. (9) 

Thus, the profiles are piecewise linear and consist of con- 
stant segments with equal densities in the two lanes and 
segments of slope O (—0) in lane A(B) with comple- 
mentary densities. Switching off the interchain particle 
exchange (O = 0), we get two identical TASEPs, which 
have, apart from the coexistence line a = [3 < 1/2, con- 
stant density profiles in the bulk. In the high-density 
phase ((3 < min{a, 1/2}), the density, being 1 — (3, is con- 
trolled by the exit rate and the profile is discontinuous at 
x = 0. In the low-density phase (a < min{/3, 1/2}), the 
density is a and a discontinuity appears at x = 1 Q, [g] . 
In the maximum current phase (a, (3 > 1/2), as we have 
already mentioned, the bulk density is 1/2 and boundary 
layers appear at both ends. On the other hand, the effect 



of symmetric lane change processes is to diminish the dif- 
ference between the local densities in the two lanes. Since 
the densities are already equal without the interaction, 
this situation is obviously not altered when switching on 
the vertical hopping processes. Consequently, the density 
profiles in the bulk are identical to that of the TASEP in 
these phases. 

This is, however not the case at the coexistence line 
a = /3 < 1/2. In the TASEP, a sharp domain wall 
emerges here in the density profile, which separates a 
low- and a high-density phase with constant densities far 
from the domain wall a and 1 — a, respectively. The 
stochastic motion of the domain wall is described b y an 
unbiased random walk with reflective boundaries [8|, [29j] , 
such that the average stationary density profile connects 
linearly the boundary densities a and 1 — a. Return- 
ing to our model, we consider first the closed system, i.e. 
a = P = 0. The profiles which fulfill the requirement 
about the continuity of the currents are depicted in Fig. 
[5] for various global particle densities g = hin^^oo 
where N is the number of particles in the system. Here, 
the density profiles consist of three segments in general. 
In the middle part of the system an equal-density seg- 
ment is found (Fig. [5^,b,d). This region is connected 
with the boundaries by complementary-density segments 
on its left-hand side and on its right-hand side, which 
are continuous at x = and x = 1, respectively. In both 
lanes, the density profile is continuous at one end of the 
equal-density segment and a shock is located at the other 
one, such that the two shocks are at opposite ends. The 
density in the equal-density region (and at the same time 
the location of the shocks) depend on the global particle 
density. At g — 1/2, the equal-density segment is lacking 
if < 1 (Fig. [2J:) and the profiles are linear if = 1. 




x 10 x 1 



FIG. 2: Density profiles of the closed system (a = /3 = 0) 
with £1 = 0.5 for different global densities: (a) g — 0.26, (b) 
g = 0.44, (c) g = 0.5 and (d) g = 0.74. The thin solid and 
dashed lines represent the flow field of the differential equation 
(J5]) corresponding to the complementary-density and equal- 
density solutions, respectively. The thick solid and dashed 
lines are the density profiles p(x) and n(x), respectively. 
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If particles are allowed to enter and exit from the sys- 
tem at the boundaries, i.e. < a = (3 < 1/2, the total 
number of particles is no longer conserved. Nevertheless, 
we expect that the stationary density profiles averaged 
over configurations with a fixed global density g, p g (x) 
and TT g (x), can still be constructed from the solutions ([8|) 
and ([5]) of the hydrodynamic equations. These profiles 
are similar to those of the closed system and the only dif- 
ference is that the complementary-density segments fit to 
the altered boundary conditions p g (0) = 7r e (l) = a and 
p e (l) = 7r e (0) = 1 — a. This is indeed the case in the limit 
a — > 0. Here, the injections and removals of particles at 
the boundary sites, which modify the global density, are 
infinitely rare, such that the system has always enough 
time to relax, i.e. to adjust the density profiles to the 
slightly altered global density. As long as the shocks 
are not in the vicinity of the boundaries, the densities 
at the boundary sites are independent from the global 
density, which influences only the position of the equal 
density segment. Therefore, the stochastic variation of 
the global density g(t) is described by a homogeneous, 
symmetric random walk in the interval [0, 1] with reflec- 
tive boundary conditions. For finite a, we can give only a 
heuristic argument why we expect that the fluctuations 
of the global density are quasistationary in the above 
sense. In the stationary state, the center of the mass of a 
small instantaneous local perturbation propagates with 
a velocity v(x) = 1 — 2p(x) j29|, |43|, which changes sign 
at p = 1/2. In the complementary-density segments, 
the perturbations in the density, which come from the 
fluctuations of the boundary reservoirs, are thus driven 
toward the equal-density segment with a finite velocity. 
The characteristic traveling time of the perturbation, as 
well as the time scale related to the lane change processes 
in a finite system of size L is O(L). The relaxation time of 
the perturbation is thus expected to be 0{L). However, 
the random walk dynamics of the global density implies 
that the time scale of a finite change in the global den- 
sity is 0(L 2 ), which is large compared to the relaxation 
time, thus the density profile has enough time to follow 
the instantaneous global density. The fluctuating global 
density g(t) is thus expected to be a symmetric random 
walk with reflective boundaries at a and 1 — a. In the 
stationary state, the global density is therefore homoge- 
neously distributed in the interval [a, 1 — a]. 

On the other hand, one can easily calculate that if the 
position of the shock in lane A is x s , the global density of 
particles in the system is g(x s ) = (1 — x s )[A(x s ) + f2] +a 
for x s > 1/2, where we have introduced the (position 
dependent) height of the shock: A(x s ) = 2£l(x s — 1) + 
1 — 2a. Note that, as opposed to the single lane TASEP, 
this relation is no longer linear, therefore the probability 
distribution of the position of the shock is not uniform 
in the steady state. 

With these prerequisites, the steady-state density pro- 
file in lane A can be easily calculated by averaging p e (x) 
over the steady-state distribution of the global density: 

p( x ) — i -2a la " Pg( x )dg- Skipping the details of the 
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FIG. 3: Density profiles in lane A at the coexistence line at 
a = (3 — 0.1, obtained by numerical simulation for different 
system sizes and for different values of fl: a) Q — 0.4, b) 
fl = 0.8, and c) Q — 1.2. In the case Q = 0.8, the height of the 
shock at x s = 1/2 is zero. The solid curves are the analytical 
predictions in the hydrodynamic limit. The thick solid lines 
represent the profile p e (x) at global density q — 1/2. 



straightforward calculations, we shall give the profile p(x) 
in the interval k < x < 1, whereas for < x < i, it is 
obtained by the help of the relation p{x) = 1 — p(l — x), 
that follows from eq. © and 0. The density profile in 
lane B can be calculated by making use of eq. ([7]) , which 
implies tt{x) = 1 — p(x). The cases corresponding to the 
different signs of A(|) must be treated separately. For 
> 0, we obtain 

P[X) = 2 + l-2a A( 2 } " °' (10) 
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which is a third-degree polynomial of x. If A(^) > 0, 
the second derivative of p(x) is discontinuous at x = h. 
In the limit -> 0, the linear profile of the TASEP at 
the coexistence line is recovered. If A(i) = 0, eq. (JTOJ) 
simplifies to 



4fi(x-i) 3 



A( 5 ) = 0, 



(11) 



which is everywhere analytic. For A(-|) < 0, the profile 
is constant in the interval i < x < |A(^)|/(2f2) + |, 
where p(x) = 1/2, while it is given by eq. (| 1 Ojl in the 
interval |A(|)|/(2Q) + \ < x < 1. These curves, as well 
as results of Monte Carlo simulations for finite systems 
of size L = 64, 128 and 256 are shown in Fig[3l In the nu- 
merical simulations, after waiting a period of 10 6 Monte 
Carlo steps in order to reach the steady state, we have 
measured the local occupancies every 10 Monte Carlo 
steps during a period of 5 • 10 9 steps. For increasing L, 
the properly scaled profiles seem to tend to the analytical 
curves expected to be valid in the continuum limit. 



IV. ASYMMETRIC LANE CHANGE 

In this section, the stationary properties of the model 
are investigated in the case K 1. Due to the sym- 
metries of the system, we may restrict ourselves to the 
investigation of the part K < 1, j3 < a of the parame- 
ter space, which is related to the remaining part through 
eqs. (|6]) and ([7]). Analyzing the solutions to the hydro- 
dynamic equation ([5]), one can construct the phase di- 
agram in the four-dimensional parameter space spanned 
by a, j3, £Ia and fls- Two representative two-dimensional 
cross sections of the parameter space at fixed lane change 
rates are shown in Fig. [4] and [5l As can be seen, the 
phase boundaries are symmetric to the diagonal, which 
is a consequence of eq. (|7|). On the other hand, the 
phase diagrams are richer compared to that of the sym- 
metric model: Besides the phases where the profiles are 
continuous in the interior of the system, the asymmetry 
in the lane change kinetics leads to the appearance of 
phases where one of the lanes contains a localized shock 
in the bulk. This is reminiscent of the shock phase of the 
single lane TASEP with Langmuir kinetics. As a new fea- 
ture, the position of the shock may vary discontinuously 
with the boundary rates here, when the so-called discon- 
tinuity line is crossed. The coexistence line, where coher- 
ently moving delocalized shocks emerge in both lanes, is 
still present, however, it is shorter than in the symmetric 
case and the shocks walk only a shrunken domain. The 
subsequent part of the section is devoted to the detailed 
analysis of these findings. 



A. Density profiles 

We start the presentation of the results with the de- 
scription of the density profiles in the phases below the 
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FIG. 4: Phase diagram at Qa = 2, Q.b = 0.4. Phase bound- 
aries are indicated by solid lines. Letters L,H and S refer to 
low-density, high-density and localized shock phase, respec- 
tively; the first(second) letter refers to lane A(B). The thick 
solid line indicates the coexistence line. At the dashed lines, 
the function J(a, j3) is nonanalytic. 




FIG. 5: Phase diagram at Q,a = 2, Qb = 1. The dotted curves 
are the discontinuity lines, which terminate at the points in- 
dicated by the full circles. 



diagonal a = (3 of the two-dimensional phase diagrams. 

If the exit rate is small enough, the densities in the 
bulk exceed the value 1/2 in both lanes (Fig. [5^). Both 
profiles are continuous in the bulk and at the exits, i.e. 
linLj—,! p{x) = lim^^o 7t(.t) = 1 — 6, but they are dis- 
continuous at the entrances, i.e. Urn x _,.Q p(x) ^ a and 
lim x ->x ir(x) 7^ a, which signals the appearance of bound- 
ary layers on the microscopic scale. The profiles p(x) 
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and tt(x), as well as the current, depend exclusively on 
(3, while a influences only the boundary layers at the 
entrances. This situation is observed also in the high- 
density phase of the TASEP with Langmuir kinetics 
[iH , Hq . therefore we call this phase H-H phase, refer- 
ring to the high density in both lanes. 



a) 






c) 


1 


















FIG. 6: Density profiles in lane A (lower line) and lane B (up- 
per line) for the parameters £Ia = 2, Qb = 0.4, a = 0.45 and 
for various exit rates: a) /3 = 0.15, b) /3 = 0.25, c) /3 = 0.35 
and d) /3 = 0.45. The profiles in panels a)-c) were obtained 
by numerically solving eq. ([5]), whereas those in panel d) are 
the analytical curves in eq. (|14p . Results of numerical sim- 
ulations (dotted lines) obtained for system size L = 10000 
by averaging the occupancies in a period of 10 7 Monte Carlo 
steps in the steady state are hardly distinguishable from these 



In the phase denoted by S-H in the phase diagram, p(x) 
and ir(x) are continuous at the exits, similarly to the fi- 
ll phase, however, the discontinuity in p(x) is no longer 
at x = but it is shifted to the interior of the system 
(Fig. Wp)- Thus ImXs-^ p{x) = a holds and a shock is 
located in lane A at some x s (0 < x s < 1). Therefore 
this phase is termed S-H phase, where the letter S refers 
to the shock in lane A and letter H refers to the high 
density (tt(x) > 1/2) in lane B. The function n(x) is 
discontinuous at x = 1, i.e. Um^^i tt(x) ^ a and it is not 
differentiable (although continuous) at x s . Since both 
p(x) and 7r(x) are continuous at x = 0, the total current 
is given by 



J = a(l - a) - 6(1 - b) 



(12) 



in this phase. The profile p(x) at fixed a and (3 can 
be computed by substituting the current calculated from 
eq. (fT2"|) into the differential equation ([5]). Then, the 
solutions propagating from the left-hand and the right- 
hand boundary, i.e. the solutions pi(x) and p r (x) fulfill- 
ing the boundary conditions pi(0) — a and p r (l) = 1 — b, 
respectively, are calculated numerically. Finally, the 
position of the shock x s is obtained from the relation 
Pi{x s ) = 1 — p r (x s ), which is implied by the continuity of 



the current in lane A. Once p(x) is at our disposal, 7r(a;) 
can be calculated from eq. (g]). 

Apart from the discontinuity line to be discussed in the 
next section, the position of the shock x s varies continu- 
ously with the boundary rates in the S-H phase. Fixing a 
and reducing (3, x s is decreasing and at a certain value of 
f3, (3 = (3h {en), the shock reaches the left-hand boundary 
at x = 0. At this point, the right-hand solution p r (x) 
extends entirely to the left-hand boundary and a further 
increase in (3 drives the system to the H-H phase. The 
phase boundary (3h{o) between the S-H and the H-H 
phase is thus determined from the condition x s = or, 
equivalently, p r (l) = 1 — a. When (3 is increased along 
a vertical path in the phase diagram at a fixed a, x s in- 
creases and for a < p\, where the constant p\ will be 
determined later, the path hits the coexistence line be- 
fore the shock would reach the right-hand boundary (see 
Sec. III. D). Increasing [3 along a path at some a > p±, 
the shock reaches the right-hand boundary at x = 1 for a 
certain value of f3, (3 = /3l(«) and the path leaves the S- 
H phase. At the phase boundary, the left-hand solution 
pi{x) extends to the whole system and pi(l) = b must 
hold. 

Crossing the phase boundary (3l{o), the L-H phase is 
entered, where letter L refers to the low density in lane 
A since here, p(x) < 1/2 and ir(x) > 1/2 hold in the bulk 
(Fig. [BJ;). In this phase, p{x) and ir{x) are discontinuous 
at x — 1, whereas they are continuous at x = hence the 
current is given by eq. (|12p . In the part of the L-H phase 
where the current is zero, i.e. if a = [3 or a, (3 > 1/2, the 
profiles can be calculated analytically. The equal-density 
solutions of eq. (O are 



ln[p e (a;)(l 



l)x + const, 



(13) 



Pe (x))] = n A {K 

ir e (x) = p e (x), 
whereas the complementary-density solutions read as 

— In \pc(x) — pi I In \pc(x) — P2I = 2QaVKx + const, 

P2 Pi 

1 



TT c (x) 



Pc(x), (14) 



where the constants p\ = 1/(1 + K~ 1 / 2 ) and P2 = 1/(1 — 
A' -1 / 2 ) are the roots of the equation Sa(p, 1 — p) = 0. 
There is, furthermore, a special complementary-density 
solution with constant densities: 



Pc(x) 
ir c (x) 



Pi, 

I- pi. 



(15) 



In the part of the L-H phase where J = 0, the profiles 
are given by the complementary-density solution which 
fulfills the boundary conditions p c (0) = a and 7r c (0) = 
1 - a (Fig. HH). 



B. Phase boundaries and the discontinuity line 

In the S-H phase (and the L-H phase), the profiles p(x) 
and 7r(a;), as well as the current are independent of a if 
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a > 1/2. Here, lim^^o p{%) = 1/2 and a influences only 
the microscopic boundary layer, as we argued in Sec. II. 
As a consequence, the phase boundaries (3h{ol) and 
are horizontal lines in the domain a > 1/2 (see Fig. 2] 
and [5]) and we may restrict the investigation of the phase 
boundaries to the domain a < 1/2. 

Although we cannot give an analytical expression for 
the density profiles in general, some information can be 
gained on the phase boundaries of the S-H phase by in- 
vestigating the constant solutions of the hydrodynamic 
equations. A constant solution p(x) = r, tt(x) = p 
must obey Sa(t,p) = 0, otherwise the spatial deriva- 
tives d x p(x) and d x ir(x) would not vanish in eq. ((2]). 
On the other hand, the constants satisfy the equation 
J = r(l — r) — p(l — p), where J is determined by the 
boundary rates via eq. (fT2|) . Eliminating p yields that r 
is given by the implicit equation: 



-(r - 1) 



1 



K 



(K + (1- K)rf 



= J(a,f3). (16) 



In the S-H phase, this equation has two roots 
rx(J(a, (3), K) and rz{J{a,[3),K) (shortly r% and r 2 ) in 
the interval [0, 1]. One can check that for the larger root 
r 2 , the relation 1/2 < 1 — (3 < r 2 holds, whereas the 
smaller one, r\, may be larger or smaller than 1/2. 

First, we examine the phase boundary separating the 
L-H phase and the S-H phase. One can check that at this 
boundary line, r\ < 1/2 holds. Moreover, it follows from 
eq. © that > if < p{x) < n, and ^1 < 

if t*i < p(x) < 1/2. Thus, the line p(x) — r*i behaves as 
an attractor for the solutions pi(x) propagating from the 
left-hand boundary x = if pi(0) = a < 1/2, meaning 
that pi (x) approaches monotonously to r\ as x increases 
and hm x _ >00 (a;) = r\. Since pi{x) is monotonous and 
pi(0) — a, as well as pi{l) — (3 hold at the phase 
boundary, at the common point of the boundary line and 
the diagonal a = /3, pi{x) must be a constant function 
pi{x) = a. This, however, implies that a must coincide 
with r%. The endpoint of the boundary line (3l{o>) is 
therefore at a = j3 = n(J = 0,K) = 1/(1 + if~ 1/2 ), 
which depends only on K . 

As opposed to this point, the whole function /3i(a) 
depends both on VIa and fig. Nevertheless, we can 
find an analytical expression for /3i(a) in the limit K = 
const, VLa — > oo. We can see from eq. ([5]) that the 
derivative dp ^ is proportional to flA for a fixed K. 
As a consequence, the larger D,a is the more rapidly 
pi{x) tends to r\. Thus, in the limit specified above, 
limn^^o^/^l) — ri) = 0, which leads to — r\. Sub- 
stituting this into eq. (|16p . we obtain for the inverse 
of the boundary curve Ploo{oi) in the limit K = const, 
SIa — * oo: 



\ 



1 - 4/?(l - (3) 



K 



[K+{1-K)0T 



This curve is plotted in Fig. [71 The phase boundaries 
obtained by integrating eq. ([5]) numerically for finite lane 
change rates tend rapidly to this limiting curve. 




.(17) 



a 



FIG. 7: Phase boundaries of the S-H phase obtained by inte- 
grating eq. <(Sj numerically for K — 1/5 and for various Qa- 
The solid lines are the limiting curves: /3_Loo(a), given by eq. 
(fT7|) and 0Hoo(ce), given solely by eq. ([18} since K < K* . 

Next, we turn to examine the boundary curve between 
the S-H and the H-H phase, /?#(a). Along this line, 
p r (0) = 1 — a and p r (l) = 1 — (3 hold, and the roots of 
eq. ([in]) are arranged asri < 1— a < l — (3<r2- One can 
show that the line p(x) — r 2 is an attractor for the solu- 
tions pi{x) which start at x = if pi(0) > max{l/2,ri}. 
Moreover, if r\ > 1/2, the line p(x) = r\ repels the so- 
lutions pi{x) starting from x = if r\ < pi{0) < r 2 , 
or, in other words, the solutions p r {x) propagating from 
the right-hand boundary, for which n < Pr(l) < T2, are 
attracted by the line p(x) — r\. When the diagonal is 
approached along the phase boundary /^(a), the pro- 
file p r (x), being monotonous, must tend to the constant 
function p r {x) = 1 — a, as well as tt(x) since the current 
is zero at the diagonal. However, equal densities in the 
two lanes are possible for K ^ 1 only if the density is one 
(or zero) , therefore the boundary line must approach the 
diagonal at a = 0. This is in accordance with the fact 
that r 2 = 1 if J = 0. Thus, we obtain lim Q _>o (3l{o) — 0, 
independently of the lane change rates. 

Similarly to /3z,(a), the location of the whole boundary 
line /3h(c() depends on both lane change rates (see Fig. [7] 
and[H]) and we can give an analytical expression for (3h (a) 
again only in the limit K = const, Ha — > 00. As afore- 
said, the profile given by p r (x) lies between two attrac- 
tors, p{x) = r\ and p(x) — r 2 , to which the solutions tend 
in the limit x — > —00 and x — * 00, respectively. When £Ia 
is increased (such that K is fixed), then /3i(a) at a fixed 
a decreases. Thus, the current is increasing and the two 
roots of eq. (|16|) . r± and r 2 are coming closer. Keeping 
in mind that n < 1 — a < 1 — /3 < r 2 , there are now two 
possible cases. Depending on the value of a, either the 
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FIG. 8: Phase boundaries between the S-H phase and the 
H-H phase obtained by numerical integration of eq. (J5]l for 
K = 1/2 > K* and for several values of £Ia- The solid lines 
are the limiting curves given in eq. (|18[l and eq. (I19p . 



gap between 1 — and r 2 or the gap between rq and 1 — a 
vanishes first. In other words, in the former case, the pro- 
file is attracted to the line p(x) — r 2 at x = 1 in the limit 
Qa — > oo, i.e. limo A _, 00 (p r (l) — r 2 ) =0, whereas in the 
latter case it is attracted to the line p{x) = r± at x = 
(provided that n > 1/2) , i.e. limn A _ >oo (/j r (0) — ri) = 0. 
In the first case, substituting r = 1 — (3 into eq. (|16[). 
we obtain for the inverse of the limiting curve b Hoo (a) 
in terms of 0Loo(oe) 



[p b Hoo ]- 1 (P) = [pLoo}- 1 (l-(3), 

whereas in the second case, r — 1 — a yields 

/3&oo(«) = 



(18) 




4a(l - a)if 



[if +(l-if)(l-a)r 



(19) 



These curves are plotted in Fig. [5] The phase boundary 
in the limit if = const, Qa — > oo is given by 



/feoo(a) = max{0% oo (a),0 b Hoo (a)}. 



(20) 



The value a* at which the functions 0%(ct) and h H {a 
intersect varies with if . If if 
0.21132 . . . , while a* = 1/2 if if is equal to 



1, a* tends to = 

2V3 



if* = 1 + V2 - a/2(1 + V2) = 0.21684 . 



(21) 



Thus, for if < K*, the limiting curve of Ph(<x) is given 
by eq. (HHJ) alone, otherwise it is composed of eq. l[TB]) 
and eq. (fT9|) as given by eq. (|20)l . 

Although the curve Phoo(&) gives the phase boundary 
line only in the limit K = const, £Ia —> oo, we show that 
for K > K* and for large enough Qa, S!a > Cl A (K), the 



phase transition point at a = 1/2 is given exactly by eq. 
d, i-e- %(l/2) = P a Hoo (l/2) - j^. 

In order to see this, we discuss first the possible ap- 
pearance of a discontinuity line in the S-H phase, at 
which the position of the shock x s changes discontinu- 
ously. If r\ — 1/2, one can see from eq. © that for 

the spatial derivative of the profile, \\m p ^i/ 2 dp }^ > 
holds and the left-hand and right-hand solutions pi(x) 
and p r {x) may propagate as far as the line p(x) = 1/2. 
If pi(x\) — 1/2 and p r (x2) = 1/2 hold for some x% and 
X2, such that < x\ < x% < 1, then the left-hand and 
right-hand solutions are connected by a constant seg- 
ment p(x) — 1/2 in the interval [aq,^] and the pro- 
file is continuous (Fig. [9]:). Substituting r = 1/2 into 





FIG. 9: Density profiles for parameters Qa = 2, Qb — 1 
and a) a = 0.4, /3 = 0.29, where n is slightly above 1/2; 
b) a = 0.4, f3 = 0.32, where n is slightly below 1/2; c) 
a = 0.4, P » 0.305, where n = 1/2. d) The endpoint of the 
discontinuity line at a ~ 0.297, /3 w 0.237. 

eq. (I16p we obtain the equation of the discontinuity line: 



a{l-a)-0{l-P) = { 



l-K 



- , along which the current 



is constant. Solving this equation for 0, we obtain 

0d{a) 




Mi 1 - a) +4 



1 - if 
2(1 + if) 



(22) 



This curve is shown in Fig. [5l When at an arbitrary 
point of this line, is infinitesimally decreased, then rq 
exceeds the value 1/2 and a shock appears at x\ with an 
infinitesimal height (Fig. [5^). Conversely, an infinitesi- 
mal increase in decreases rq below 1/2 and an infinites- 
imal shock appears at x 2 (Fig. Od). Thus, when this line 
is crossed, the position of the shock jumps from x\ to 
X2- At the point of the discontinuity line at a = 1/2, 
Xi = holds and an infinitesimal increase (decrease) in 
drives the system to the H-H (S-H) phase. There- 
fore the point of this curve at a — 1/2 coincides with 
the phase boundary, i.e. 0d(l/2) = /3ff(l/2). On the 
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other hand, comparing eq. (|19|) and eq. (|22)) . we see that 
(3 d (l/2) - a Hoo (l/2) = j^. The function (3 a Hoo (a) thus 
gives the phase transition point at a — 1/2 exactly, pro- 
vided the curve (3d{a) is located in the S-H phase. That 
means if K > K* and SIa > ^a(-^0; where W A (K) is 
the value of id a at which (3h (1/2) first reaches the value 
when VLa is increased from zero. For example, for 



K 
l+K 



K = 1/2, we have found 0^(1/2) » 0.719 (Fig. |gj). We 
emphasize, however, that the discontinuity line is lacking 
if K < K* or K > K* but Q A < tt* A ( K )- 

If K > K* and Ha > n A (K), then, at the transition 
point at a = 1/2, we have x\ = 0; Qa influences only 
#2 and X2 — > if — > f2^(i4T). Moving away from the 
point at a = 1/2 along the discontinuity line, the length 
of the constant segment X2—x\ is decreasing and vanishes 
at a certain point. Here, the density profile becomes 
analytical and the discontinuity line terminates (see Fig. 
|SJi and Fig. [5]). The position of this endpoint depends 
on Qa and it is moving toward smaller values of a for 
increasing fi^; in the large £Ia limit, it tends to the point 
of intersection of (3d(ct) and (3 b Hoo , whereas in the limit 
ft A — » fl A (K), the abscissa of the end point tends to 1/2. 

Finally, we mention that another special curve in the 
S-H phase is defined by the equation Sa(o>, 1 — (3) = 0. 
At this curve the left-hand solution pi(x) is constant, 
therefore both p(x) and nix) are constant in the interval 
[0,x.]. 




FIG. 10: The current as a function of /3 along the line a = 1/2 
for K = 1/5 (a) and K = 1/2 (b). The solid line is the current 
in the limit Qa — > oo, the other curves from bottom to top 
correspond to £Ia = 0.25, 0.5, 1, 2, 4, respectively. 



in the domain K < K* . For K < K* , the maximal cur- 
rent is thus J msx (K) = 1/4-/3^ (1/2) (1 -/W(l/2)). 
Since /3jjoo(l/2) decreases monotonously with decreasing 
if, the current is maximal at X = 0, where /3ffoo(l/2) = 
2 ~ 4 ^ and J m ax(0) = 1/8. The value 1/8 is thus an up- 
per bound for the total current and J — > 1/8 in the limit 
fl A -► oo if a > 1/2, (3 = and VL B = 0. 



C. Current 

We have seen that, apart from the H-H (and the L-L 
phase), the current is given by eq. (fl"2")) . It is zero at 
the line a = f3 and in the domain a, (3 > 1/2 and it 
is independent of a (f3) in the H-H (L-L) phase. It is 
a continuous function of the boundary rates everywhere, 
although nonanalytic at the boundaries of the H-H phase 
and the L-L phase, as well as at the lines a = 1/2 and 
(3 = 1/2 outside the H-H phase and the L-L phase, where 
the effective boundary rates a and b saturate at 1/2. 

In the following, we concentrate on the current in the 
H-H phase. Since it is continuous at the boundary line 
Ph(ch), it can be expressed in the H-H phase in terms of 
p H (a) as J{(3) = f3 H l {[3){\ - (3 H l {(3)) - (3(1 - (3). Ac- 
cording to numerical results (see Fig. [TO)) , the current 
J((3) has a maximum in the H-H phase, and for large 
Ha, the location of the maximum tends to (3 — /3jy (1/2) 
for K < K* , while for K > K*, it tends to (3* where the 
curves ^^{a) and f3 b Hoo (a) intersect. Contrary to this, 
the current is a monotonously decreasing function of (3 
in the S-H phase. We now turn to the question, at which 
parameter combination the total current is maximal in 
the steady state. For fixed K, the maximal current is 
realized in the limit VLa — > oo at (3hoo(1/2) for K < K* 
and at (3* for K > K* . Since the current is growing 
faster in the S-H phase for decreasing (3 than in the H-H 
phase above (3* (if K > K*), we conclude that the pa- 
rameter combination that maximizes the current is found 



0.02 




P 

FIG. 11: The current as a function of (3 in the H-H phase ob- 
tained by numerical integration for K — 1/2 and for different 
values of Oa- The analytical approximation given in eq. (123[) 
is indicated by solid lines. 



We close this section with the discussion of the current 
in the H-H phase when the lane change rates are small, 
i.e. £Ia,Qb "C 1. If the vertical hopping processes are 
switched off, the current is zero, hence we expect that 
for small lane change rates the current is small, as well. 
Assuming that J <§C (| — (3) 2 and expanding the right- 
hand side of eq. ^ in a Taylor series up to first order in 
J, then solving the resulting differential equation yields 
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finally 



•/ = ^—^/3(1- /?)(!- 2(3) x 



l + K' 



'1 + 2(1 + 10 



1 - 2(3 



(23) 



The details of the calculation are presented in Appendix 
A. This expression is compared to the current calculated 
by integrating eq. ([5]) numerically in Fig. 1111 Expanding 
this expression for small £l Aj we obtain 



j = 0(i-0)(i-K)nj 



i + K n 



A 



1 - 2/3 



a)- 



The current is thus in leading order proportional to Qa- 
Examining the higher-order terms in the series expansion 
of the right-hand side of eq. (JSJl , one can show that for 
arbitrary Q A , the current vanishes as J ~ 1 — K when 
K -»■ 1. 




FIG. 12: The flow field of the differential equation (0 for 
Q.A = 2, Qb = 0.5 and J = 0. The solid curves represent the 
complementary-density solutions given in eq. (|14|l . whereas 
dashed curves represent the equal-density solutions given in 
eq. (H3J. 



D. Coexistence line 

Let us assume that a given point of the section a = 
(3 < pi is approached along a path in the S-H phase. In 
this case, the position of the shock in lane A, x s , tends 
to some Xmin = £min(a, Q A , Vtg), which is somewhere 
in the bulk, i.e. < x m i n < 1. When the same point 
is approached from the L-S phase, the position of the 
shock in lane B tends to the same x m j n according to 
eq. ([7]). Thus, when the boundary line a — (3 < p\ is 
passed from the S-H phase, the shock in lane A jumps 
from a; m i n to the right-hand boundary at x = 1, where 
a discontinuity appears, while the discontinuity at x = 1 
in lane B, which can be regarded as a shock localized 
there, jumps to x m i n . So, the density profile changes 
discontinuously. Strictly at a = (3, the shocks in both 
lanes are delocalized and perform a stochastic motion in 
the domain [x m i n ,l], similarly to the symmetric model 
with K = 1 at the line a = < 1/2. 

Now, this phenomenon is investigated in detail in the 
case of asymmetric lane change. Since the current is zero 
if a — [3, the solutions of the hydrodynamic equations are 
those given in eq. (JT3J) and eq. {SJ (see Fig. Q2J). The 
argumentation about the quasistationarity of the fluctu- 
ations of the global density presented in the case K = 1 
apply to the case K ^ 1, as well. Thus, in the open 
system, the density profiles averaged over configurations 
with a fixed global density, p g {x) and %g{x), can be con- 
structed from the solutions (jT3)) and (fT4]) . 

The structure of these profiles is identical to those 
obtained in the case K = 1 (see Fig. [T3J) . Generally, 
they consist of three segments (see Fig. Q3b): An equal- 
density segment is located in the middle part of the sys- 
tem in the domain [xo, xo]. This region is connected with 
the left-hand boundary by a complementary-density seg- 
ment, which is continuous at x = 0, i.e. p c ,i(0) = ol, 



a) b) 























X mln 

c) 


d) 




FIG. 13: Density profiles in lane A (thick solid lines) and in 
lane B (thick dashed lines) which belong to a fixed global 
particle density q, at four different values of q. The flow field 
of eq. ([5)l is indicated by thin lines. The rates are Qa = 0.5, 
Qb = 0.25 and a = (3 = 0.1. 



"0/(0) = 1 — a, and with the right-hand boundary by 
another complementary-density segment, which is con- 
tinuous at x = 1, i.e. p c ,r(l) = 1 — ot, 7T C)r (l) = a. 
Each lane contains a shock, which are at the opposite 
ends of the equal-density segment (one at xo, the other 
one at xq). The location of the equal-density segment, 
as well as Xo and xo are determined by the actual global 
density. If xo = xo = %, the equal-density segment is 
lacking and p c ,i(%) and p c ,r(x) are directly connected by 
a shock at x (Fig. [T3"bi. Since xo and xo are not in- 
dependent, the shocks move in a synchronized way, and 
their motion is confined to the range [x m j n ,l]. If one 
of the shocks is at x = 1, the other one is at x m ; n , 
thus, the lower bound x m ; n is determined by the equa- 
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tion p c ,/(x min ) = 1 - p e (x min ), where p e (x) is the equal- 
density solution which fulfills p e (l) = 1 — a (see Fig. 
[13k) . The lower bound x m i n is an increasing function of 
a (see Fig|T4"|). In the limit a — * 0, x m \ n tends to zero and 
if a — ► pi, x m i n tends to 1, thus, at a = pi, the shock 
becomes localized at x = 1 and the system enters the 
L-H phase. At this point, the density profile is given by 
the special complementary-density solution: p c (x) = a, 
tt c (x) = 1 - a. 




0.5 1 



x 

FIG. 14: The shock in lane A at the left most possible position 
for Qa = 0.5, f2s = 0.25 and for different rates a(= 0). 

Similarly to the case K = 1, the stochastic variation 
of the global density g(t) is described by a bounded sym- 
metric random walk. The stationary density profile p(x) 
can be obtained from the profile p g {x) at a fixed global 
density by averaging it over g. For K ^ 1, we could not 
carry out the averaging analytically, nevertheless, we can 
gain some information on the stationary density profiles 
at x without the knowledge of the entire profiles. At the 
point x 7 the relation p g (x) = Tr e (x) holds for all g. This, 
together with the relation p(x) = 1 — n{x) following from 
eq. ([7]) implies that p{x) = -k(x) = 1/2. As it is shown 
in Appendix B, the ratio of the first derivatives of the 
stationary density profile in lane A on the two sides of 
the point x is 

dp(x-) dp(x+) = K+(l- K)p 

dx 1 dx 1 - (1 - K)p ' 1 J 

where p = Pc, a {x). In the case K < 1, po < pi < 1/2 
always holds, hence this ratio is smaller than 1 and the 
first derivative of the density profile is discontinuous at x. 
Furthermore, it is clear that the stationary density profile 
is identical to the complementary-density solution in the 
interval [0,x m i n ], since this domain is forbidden for the 
shocks. We have performed numerical simulations for 
finite systems of size L = 64, 128 and 256 and measured 
the density profiles in the same way as in the symmetric 
case at the coexistence line. Results are shown in Fig. 

ma 



1 







L=64 




l_=128 




L=256 - 





(i-0.5)/L 1 



FIG. 15: Density profiles in lane A obtained by numerical 
simulations with rates a — f3 = 0.1, Qa = 0.5, Qb = 0.25 and 
for different system sizes (thin lines) . The thick line represents 
the density profile p g (x) at global density g = 1/2 in the 
continuum limit. 

V. DISCUSSION 

In this work, a two-lane exclusion process was stud- 
ied, where the particles are conserved in the bulk and 
each lane can be thought of as a totally asymmetric sim- 
ple exclusion process with nonconserving kinetics in the 
bulk. As a consequence, the model unifies the features 
of the particle conserving and bulk nonconserving exclu- 
sion processes, as far as the dynamics of the shock is 
concerned. Namely, it exhibits phases with a localized 
shock in one lane, while the other one acts as a nonho- 
mogeneous bulk reservoir, the position-dependent den- 
sity of which is determined by the dynamics itself in a 
self-organized manner. On the other hand, the model 
undergoes a discontinuous phase transition at the coex- 
istence line, where delocalized shocks form in both lanes 
and move in a synchronized way. Here, the global density 
of particles behaves as an unbiased random walk, simi- 
larly to the TASEP at the coexistence line, however, the 
density profiles in the coexisting phases are not constant 
here. 

Although we considered throughout this work ex- 
change rates proportional to 1/ L, one may imagine other 
types of scaling. For the TASEP with Langmuir kinet- 
ics, shock localization is observed at the coexistence line 
when the creation and annihilation rates vanish propor- 
tionally to l/L a with 1 < a < 2 [20]. It might be worth 
examining whether the synchronization of shocks in the 
present model persists for lane change rates vanishing 
faster than 1/L. 

In the limit of large lane change rates, K = const, 
Ha — > co, the density profiles and the boundaries of the 
phases exhibiting a localized shock are related to the ze- 
ros of the source term in the hydrodynamic equation, 
which is generally valid for systems with weak bulk non- 
conserving kinetics 0, H^, [Z^]. However, different be- 
havior is observed in general when the lane change rates 
loa and ujb are finite in the limit L — > oo. In this case, the 
particle current from one lane to the other may be finite 
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at the boundary layers. These currents add to the incom- 
ing and/or outgoing currents at the boundaries, which 
may lead to nontrivial bulk densities even in the case 
K = 1, where the profiles are constant. 

Possible extensions of the present model are obtained 
when different exit and entrance rates or different hop 
rates in the two lanes are taken into account. Neverthe- 
less, these generalized versions are more difficult to treat 
because of the reduced symmetry compared to that of 
the present model. 



APPENDIX A 

We derive here an approximative expression for the 
current in the H-H phase in the limit of small lane change 
rates. Assuming that J -C (h — /3) 2 , we may expand the 
right-hand side of the differential equation © in a Taylor 
series up to first order in J. Integrating the differential 
equation obtained in this way yields 



F(p) 



1 



K - 1 
1 

0~K) 



ln(p - p 2 ) + J 



K 



In- 



{l-KY 
1 



In- 



P 



1 



P 



1 



1-p p 
Qax + const, 



(Al) 



APPENDIX B 

As we have seen, the positions of the shocks in lane 
A and lane B are not independent, thus, we may define 
thereby a function xo(xo), which is given implicitly by the 
equation p e (xo) = 1 — Pc,r(xo), where p e (x) is the equal- 
density solution which satisfies the condition p c i(xo) — 
p e (xa). In the following, we shall denote the density in 
lane A at xo, if the shock is located at x s , by p Xa (xo)- 
First, we notice that at the reference point Xq (xo < x), 
Px a {xo)=Pc,i{xo) holds whenever the shock in lane A re- 
sides between Xq and Xo(xo), i.e. xq < x s < Xo(xq). Sim- 
ilarly, p Xa (x {xo))=p C:r {xo) if x Q < x s < x (x ). When 
the shock is outside this interval (x a ^ [xo,£o])j then 
Pxs i x o)—Pe(xo), where p e (x) is some equal-density solu- 
tion determined by the global density and p Xa (xo) > 1/2 
or p X3 (xo) < 1/2 if x s < xo or x s > xo(xo), respectively. 
As a consequence of the particle-hole symmetry, the re- 
lation px.(xo) = 1 - p Xe ( Xt )(x ) holds if x s £ [x ,Xo]. 
Furthermore, in the stationary state, the probability that 
x s < xq is equal to the probability that x s > Xo(xq) for 
any Xq < x. Therefore the contribution to the average 
profile is 1/2 when x s </ [xo,Xq]. We can thus write for 
the average densities at Xq < x and Xq(xq) > x 



p(x ) = p(x )p c ,i(x ) + [1 -p(xo)] 



where the first term on the left-hand side is just the equal- 
density solution for J = 0. The solution which obeys the 
boundary condition p(l) = 1 — j3 is F(p) — F(l — j3) = 
£Ia(x — 1). From this equation, we get for the density at 
the left-hand boundary, po = lim^^o p(^), the implicit 
equation: 



F(p ) - F(l - /3) = -n A . 



(A2) 



On the other hand, we have another relation between J 
and po- 



p(x ) = p(xo(x Q ))p^ r (x Q ) + [1 - p(x (x ))]-, (Bl) 

respectively, where p(xo) is the probability that the shock 
in lane A resides in the interval [xo, xo(xo)], and Xo(xo) is 
the inverse function of xo(a;o)- For the spatial derivatives 
of the densities at xo and xo(xo), we get 

p'(x ) =p'(x ) (pc,i( x o)- ^ +p(xo)p'cAxo), 



J = po(l-Po)-P(l-P), 



(A3) 



p'(x ) =p'(x )x (x ) (pc,i(xo) - r ) +p{xo(xo))p'(x ), 



(B2) 



thus, we have closed equations for current. Note that 
the leading order term on the left-hand side of eq. (|A2|) . 
which comes from the difference of the leading terms of 
F(p) evaluated at po and 1 — j3, is O(J), while the next- 
to-leading contribution in eq. (|A2[) is 0(J 2 ). Expressing 
Po from eq. (|A3[) and expanding it for small J, we get 

po = 1 — P — jz^g — (1-2^)3 + 0(J 3 ). Substituting this 
expression into eq. (IA2j) gives an implicit equation for J. 



where the prime denotes derivation. Using that p(x) = 
and p c ,i(x) = 1 — p c , r (x), we obtain for the ratio of the 
left- and right-hand side derivatives at x: 



p'(x-)/p'(x+) = |x (x)|. 



(B3) 



Expanding the functions p Cj /(x), p c>r (x) and p e (x) in Tay- 
lor series up to first order in x around x = x, we obtain 



Assuming that J <C j3 and expanding the terms contain- 
ing J in this equation in Taylor series, we obtain finally 



|x (x)| 



P'c,l( X ) ~ Pe( x ) 
P'e( x ) ~ P'cA x Y 



(B4) 



e 



1 + K 



e 2 + o(e 3 ) 



1-2/3' 



(A4) 



Using eq. @, these derivatives can be given in terms of 
po = p CjCt (x) and we arrive at eq. (|25[) . 



where O 



For small 0, which 



(l-2/3)/3(l-/3)(l-K)- 

amounts to Qa <C 1 — 2/3, we get a good approxima- 
tion for the current by solving this quadratic equation 
and arrive at eq. 
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